Método de la transformada inversa

Daremos por hecho que podemos generar variables aleatorias independientes e identicamente distribuidas con distribución Uniforme$(0,1)$.

Empezamos por la distribución más sencilla posible: $\text{Bernoulli}(p=0.5)$.

Teniendo acceso a una variable aleatoria Uniforme$(0,1)$; ¿cómo simularían la variable aleatoria $X\sim \text{Bernoulli}(0.5)$?

En efecto $X=\mathbb{1}_{\left\{ U \leq 0.5 \right\}}$ se distribuye $\text{Bernoulli}(p=0.5)$

En general si tenemos una distribución discreta $X$ tal que para $x_i\in \left\{ x_1,x_2,\ldots ,x_n\right\}$

\begin{align*} \mathbb{P}\left[X=x_i \right]=p_i \end{align*}

con

\begin{align*} \sum_{i=1}^n p_i =1. \end{align*}

¿Cómo se puede simular $X$ haciendo uso de una variable aleatoria $U\sim\text{Uniforme}(0,1)$?

En efecto podemos definir

\begin{align*} X=\begin{cases} x_1 & \text{ si } U\leq p_1 \\ x_2 & \text{ si } p_1 < U \leq p_1 + p_2 \\ \vdots \\ x_n & \text{ si } p_1 + p_2 + \ldots + p_{n-1}< U \end{cases} \end{align*}

La metodología anterior puede ser generalizada haciendo uso de $F^{\leftarrow}$, la inversa generalizada de la función cumulativa de distribución $F$ asociada a la variable alaeatoria $X$ que se quiere simular.

Recordamos que la función de distribución cumulativa de una variable aleatoria $X$ está dada por

\begin{align*} F(x) = \mathbb{P}\left[ X \leq x \right]. \end{align*}

Una función de distribución cumulativa debe ser no-decreciente, continua por la derecha y tal que

\begin{align*} \lim_{x\to-\infty}F(x)=0 \end{align*}

y \begin{align*} \lim_{x\to \infty}F(x)=1. \end{align*}

Definimos la inversa generalizada de una función cumulativa de distribución como

\begin{align*} F^{\leftarrow}(u)=\inf\left\{ x \, : \, F(x) \geq u \right\} =\min\left\{ x \, : \, F(x) \geq u \right\}, \quad 0<u<1. \end{align*}

Proposición

Dada $F$ una función de distribución cumulativa:

• a) $u\leq F(x) \iff F^{\leftarrow}(u) \leq x$.
• b) Si $U\sim\text{Uniforme}(0,1)\implies F^{\leftarrow}(U)$ tiene función cumulativa de distribución $F$ ( es decir $\mathcal{L}(X)=\mathcal{L}(F^{\leftarrow}(U))$).
• c) Si $F$ es continua $\implies F(X)\sim \text{Uniforme}(0,1)$.

Ejemplo

Si $X\sim \text{Exponencial}(\lambda) \implies F(x)=1-e^{-\lambda x}\implies F^{\leftarrow}(x)=F^{-1}(x)=-\frac{\log(1-x)}{\lambda}$, en pracica se usa $X=-\frac{\log(U)}{\lambda}$ para simular variables aleatorias exponenciales.

using Plots, RandomNumbers.Xorshifts
r = Xoroshiro128Plus(0x1234567890abcdef)

Xoroshiro128Plus(0xe7eb72d97b4beac6, 0x9b86d56534ba1f9e)

n = 1000
x = -log.(rand(r,n));

plot(p)

plot(p)

Simulación de variables aleatorias condicionadas a un intervalo

Si se condiciona $X$ a $X\in (a,b)$ para $a<b$, entonces para $x\leq b$

\begin{align*} \mathbb{P}[X < x | X\in(a,b)] = \frac{ \mathbb{P}[a < X < \min\left\{ x , b\right\} ] }{\mathbb{P}[a<X<b]} = \frac{ F(x) -F(a)}{F(b)-F(a)}. \end{align*}

Se sigue que $F^{-1}(F(a)+U(F(b)-F(a)))\sim X|X\in (a,b)$.

Problema con el método de la transfromación inversa

• Si $F^{\leftarrow}$ no es analíticamente conocida se tiene que usar aproximaciones a esta función que suelen ser costosas de evaluar. Aún más, la inversa generalizada de $F$ puede tener dependencia compleja de los parámetros en el sentido de que para cierto valor de estos la inversa es menos costosa que para otros.

Método de aceptación rechazo

Sea $X$ una variable aleatoria continua, con función de densidad de probabilidad $f(x)$, de la cual deseamos generar simulaciones. La idea del método de aceptación-rechaza es usar una variable aleatoria auxiliar $Y$, de la cual podemos simular y un criterio de aceptación de la variable aleatoria tal que la distribución de $Y$ condicionada al criterio de aceptación coincida con la distribución de interés asociada a $X$.

En particular consideraremos $Y$ una variable aleatoria auxiliar continua con densidad de probabilidad $g$ y una constante $0<M\in \mathbb{R}$ tales que $f(x) \leq Mg(x) \, \forall x\in \mathbb{R}$.

Algoritmo de Aceptación-Rechazo

1) Simule $Y\sim \mathcal{L}(g)$ y $U\sim \text{Uniforme}(0,1)$.

2) Acepte $X=Y$ si se cumple que $$U\leq \frac{f(Y)}{Mg(Y)};$$ en caso contrario regrese a 1).

Proposición

• a) Sea $A$ el evento de aceptar en algoritmo de aceptación-rechazo, entonces $$\mathbb{P}\left[ A \right]=\frac{1}{M}.$$
• b) Si $X$ está generado por el algoritmo de aceptación-rechazo, entonces $$\mathbb{P}[X\leq x]=\int_{-\infty}^x f(z)\mathrm{d}z=F(x).$$

Intuición usando el teorem fundamental de la simulación

Observemos que una función de densidad de probabilidad $f$ siempre se puede escribir como

\begin{align*} f(x) = \int_0^{f(x)}\mathrm{d}u \end{align*}

por lo que $f$ es la densidad marginal de la distribución conjunta

\begin{align*} (X.U) \sim \text{Uniforme}\left\{ (x,u) \, : \, 0<u<f(x) \right\} \end{align*}

Teorema fundamental de la simulación

$$ X\sim \mathcal{L}(f(x)) $$

es equivalente a $$ (X,U)\sim \text{Uniforme}\left\{ (x,u) \, : \, 0<u<f(x) \right\} $$

Simular una dupla $(X,U)$ como en el teorema fundamental es sencillo si simulamos $X\sim \mathcal{L}(f)$ y $U | X \sim \text{Uniforme}(0,f(X))$, sin embargo esto no es de utilidad cuando la simulación de interes corresponde a $\mathcal{L})(f)$. Alternativamente, una metodología de aceptación-rechazo interesante en este contexto consiste en considerar $(Y,V)\sim \text{Uniforme}(M)$ con $M$ un conjunto tal que $\left\{ (x,u) \, : \, 0<u<f(x) \right\}\subseteq M$ y aceptar si la dupla $(Y,V)\in \left\{ (x,u) \, : \, 0<u<f(x) \right\}$,

Para ilustrar tal método de aceptación-rechazo consideramos una densidad $f$ tal que $$ \int_a^b f(x)\mathrm{d}x=1 $$ para $a<b$ números reales. Entonces podemos considerar $m=\max\left\{ f(x) \, : \, a \leq x \leq b \right\} $ y $M=\left\{ (y,v) \, : \, 0<v<m \wedge a<y<b \right\}$. Usamos el siguiente algoritmo de simulación.

• 1) Simular $Y\sim \text{Uniforme}(a,b)$ y $V\sim \text{Uniforme}(0,m) $ para generar una dupla $(Y,V)$
• 2) Si $V\leq f(Y)$ aceptamos $X=Y$ en caso contrario regresar al paso 1).

plot(p)

plot(p)

Anteriormente consideramos el conjunto $M=\left\{ (y,v) \, : \, 0<v<m \wedge a<y<b \right\}$ que nos ayudo a definir variables auxiliares distribuidas uniformemente sobre $M$. Más en general podemos considerar.

$M=\left\{ (y,v) \, : \, 0<v<m(y) \right\}$. Donde sin pérdida de generalidad podemos escribir $m(y)=Mh(y)$ para $0<M<\infty$ y $h$ una densidad de probabilidad.

Observemos que con $V | Y \sim \text{Uniforme}(0,m(Y))$ la regla de aceptación "$ X=Y \iff V\leq f(Y)$" tiene sentido para un algoritmo de simulación de $f$ si $f(x)\leq m(x)\, \forall x\in \mathbb{R}$; aún más en tal caso se tiene que $\mathcal{L}(X)=\mathcal{L}(f) \iff h=g$.

De la observación anterior para $g$ densidad de probabilidad y $M$ constante positiva tales que $f(x)\leq Mg(x) \, \forall x\in \mathbb{R}$ obtendríamos el siguiente algoritmo

1) Simule $Y\sim \mathcal{L}(g)$ y $V\sim \text{Uniforme}(0,Mg(y))$.

2) Acepte $X=Y$ si se cumple que $$V\leq f(Y)$$ en caso contrario regrese a 1).

Que coincide con el algoritmo de aceptación-rechazo al considerar que $U\sim \text{Uniforme}(0,1)\iff KU\sim \text{Uniforme}(0,K).$

plot(p)

Observación

Para el algoritmo de aceptación-rechazo es necesario poder evaluar $Mg(x)$ y deseable que esta evaluación tenga un bajo costo computacional. Más aún el algoritmo puede ser modificado para basarse en la evaluación de $cf(x)$ para alguna constante $c>0$, esto permite hacer simulación cuando la constante de normalización de la densidad $f$ es inaccessible o costosa de evaluar.

Algoritmo de Aceptación-Rechazo 2

1) Simule $Y\sim \mathcal{L}(g)$ y $U\sim \text{Uniforme}(0,1)$.

2) Acepte $X=Y$ si se cumple que $$U\leq \frac{cf(Y)}{Mg(Y)};$$ en caso contrario regrese a 1).

En este caso la probabilidad de aceptar es $\frac{c}{M}$ y se mantiene que $X\sim \mathcal{L}(f)$.

Podemos considerar

\begin{align*} f(x) & \propto e^{-x^2/2}\left( sin(6x)^2 +3cos(x)^2 sin(4x)^2 + 1 \right) \\ g(x) & \propto e^{-x^2/2} \end{align*}

y $M=5$

plot(p)

p=histogram(y_vec[accept_ind], nbins=250, normalize=true, label="")
title!("Histograma de simulaciones de f con aceptación rechazo")

plot(p)

Usar una función para "ensandwichar" en aceptación-rechazo

Si podemos evaluar de manera no costosa una función $g_l$ tal que

$$ g_l(x) \leq f(x) \leq M g_m(x) $$

con $M>0$ constante positiva y $g_m$ densidad de probabilidad tales que se puede implementar el algoritmo de aceptación rechazo. Entonce tiene sentido proponer el siguiente algoritmo cuando evaluar $g_l$ es menos costoso computacionalmente que evaluar $f(x)$.

Algoritmo de Aceptación-Rechazo-Sandwich

1) Simule $Y\sim \mathcal{L}(g_m)$ y $U\sim \text{Uniforme}(0,1)$.

2) Acepte $X=Y$ si se cumple que $$U\leq \frac{g_l(Y)}{Mg_m(Y)};$$ en caso contrario proceda a 3).

3) Acepte $X=Y$ si se cumple que $$U\leq \frac{f(Y)}{Mg_m(Y)};$$ en caso contrario regrese a 1).

Densidades Log-cóncavas

Decimos que una función de de densidad de probabilidad es log-cóncava si $h(x)=log(f(x))$ es cóncava, si $f$ es suficientemente suave esto es equivalente a que $h'(x)=\frac{\mathrm{d}h}{\mathrm{d}x}$ sea estrictamente decreciente o $h''(x)=\frac{\mathrm{d}^2 h}{\mathrm{d}x^2}<0$. Para este tipo de densidades es sencillo construir funciones que "ensandwichen" a la densidad de interés.

plot(p)

plot(p)

plot(p)

Recordamos que si $h$ es diferenciable la recta tangente a la gráfica de $h$ en $x_i$ está dado por $T(x)=h(x_i)+h'(x_i)(x-x_i)$.

Dados $x_i<x_{i+1}$, las rectas tangentes a $h$ en $x_i$ y $x_{i+1}$ se intersectan en el punto dado por

\begin{align*} z_i = \frac{h(x_{i+1}) -h(x_i) -x_{i+1}h'(x_{i+1}) +x_ih'(x_i) }{h'(x_i)-h'(x_{i+1})} \end{align*}

Por lo que dados puntos $T_n = \left\{ x_1, \ldots , x_n \right\}$ definimos para $x\in(z_{i-1},z_i]$

$$ M g_u(x) = h(x_i) + h'(x_i)(x-x_i) $$

con $i\in \left\{1, \ldots, n\right\}$, $z_0$ el mínimo del soporte de $h$, en caso de que $h$ tenga soporte acotado por debajo ó $-\infty$ en caso contrario, y $z_k$ es el máximo del soporte de $h$ en caso de que $h$ tenga soporte acotado por arriba ó $\infty$ en caso contrario.

Por otro lado definimos para $x\in [x_i,x_{i+1}]$

$$ g_l(x) = \frac{(x_{j+1}-x)h(x_j)+(x-x_j)h(x_{j+1})}{x_{j+1}-x_j} $$

y si $x<x_1$ o $x>x_n$ definimos $g_l(x)=-\infty$.

Algoritmo de Aceptación Rechazo Adaptativo

1) Incialize $n$ y $T_n$ (el conjunto puntos nudo).

2) Genere $X\sim \mathcal{L}\left( g_u(x) \right)$ y $U \sim \text{Uniforme}(0,1)$.

3) Si $$U\leq \frac{ g_l(X) }{Mg_u(X)}$$ acepte $X$; en caso contrario si $$U\leq \frac{ f(X) }{Mg_u(X)} $$ acepte $X$ y actualize $n=n+1$, $T_{n+1}=T_n \cup \left\{ X \right\}$; en caso contrario regrese a 2).

@time prueb_prev = adaptive_accept_reject_vprev([-1.1,1.1],10000,h,dh);

 29.384700 seconds (1.61 G allocations: 26.664 GiB, 14.47% gc time)

@time prueb = adaptive_accept_reject([-1.1,1.1],10000,h,dh);

  0.606691 seconds (12.94 M allocations: 319.221 MiB, 12.00% gc time)

@time randn(10000);

  0.000094 seconds (6 allocations: 78.359 KiB)

plot(p)

plot(p)

plot(p)

Ejemplo

La distribución gaussiana inversa generalizada está determinada por la siguiente función de densidad de probabilidad

\begin{align*} f_{GIG}(x)= \frac{\left(a/b\right)^{p/2}}{2\mathcal{K}_p\left(\sqrt{ab}\right)} x^{p-1}e^{-\left(ax +b/x \right)/2}\mathbb{1}_{\left\{ x >0 \right\}} \end{align*}

con $a,b>0$, $p\in \mathbb{R}$ y $\mathcal{K}_p$ la función modificada de Bessel del segundo tipo.

Observemos que

\begin{align*} f_{GIG}(x) \propto x^{p-1}e^{-\left(ax +b/x \right)/2}\mathbb{1}_{\left\{ x >0 \right\}} \end{align*}

por lo que existe $c>0$ tal que \begin{align*} h(x)=\log \left( c f_{GIG}(x) \right) =(p-1)\log(x)-\frac{ax}{2} - \frac{b}{2x} \end{align*} satisface \begin{align*} h''(x)=-\frac{(p-1)}{x^2} -\frac{b}{x^3} \end{align*} por lo que $h$ es cóncava para $p\geq1$ con $x$ restringida a $(0,\infty)$.

plot(p)

plot(p)

Ejemplo

Regresión Poisson

Supongamos que tenemos muestras $(Y_1,X_1), \ldots , (Y_n,X_n)$ tales que

\begin{align*} Y_i | a,b,X_i & \sim \text{Poisson}\left( e^{aX_i + b} \right) \\ a & \sim \text{Normal}(0,\sigma^2) \\ b & \sim \text{Normal}(0,\tau^2) \end{align*}

Si consideramos observaciones $Y_i=y_i$, la distribución posterior de $a,b|\pmb{Y}=\pmb{y},\pmb{X}=\pmb{x}$ entonces tiene densidad dada por

\begin{align*} \pi(a,b|\pmb{y},\pmb{x}) &\propto \pi(\pmb{y}|\pmb{x},a,b) \pi(a) \pi(b) \\ & \propto exp\left( a\sum_{i=1}^n y_i + b\sum_{i=1}^n y_i x_i - e^{a}\sum_{i=1}^n e^{x_i b} -\frac{a^2}{2\sigma^2} -\frac{b^2}{2\tau^2} \right) \end{align*}

Entonces

\begin{align*} log\left( \pi(a|\pmb{y},\pmb{x},b) \right) & = a\sum_{i=1}^n y_i - e^{a}\sum_{i=1}^n e^{x_i b} -\frac{a^2}{2\sigma^2} \\ log\left( \pi(b|\pmb{y},\pmb{x},a) \right) & = b\sum_{i=1}^n y_i x_i - e^{a}\sum_{i=1}^n e^{x_i b} -\frac{b^2}{2\tau^2} \end{align*}

Podemos considerar \begin{align*} h_1(a) = & a\sum_{i=1}^n y_i - e^{a}\sum_{i=1}^n e^{x_i b} -\frac{a^2}{2\sigma^2} \\ h_2(b)= & b\sum_{i=1}^n y_i x_i - e^{a}\sum_{i=1}^n e^{x_i b} -\frac{b^2}{2\tau^2} \end{align*} con lo que \begin{align*} h_1'(a) & = \sum_{i=1}^n y_i - e^{a}\sum_{i=1}^n e^{x_i b} -\frac{a}{\sigma^2} \\ h_2'(b) & = \sum_{i=1}^n y_i x_i - e^{a}\sum_{i=1}^n x_i e^{x_i b} -\frac{b}{\tau^2} \end{align*}

using RDatasets
prussian = dataset("pscl", "prussian")

┌ Info: Recompiling stale cache file /home/workstation/.julia/compiled/v1.2/RDatasets/JyIbx.ji for RDatasets [ce6b1742-4840-55fa-b093-852dadbb1d8b]
└ @ Base loading.jl:1240

280 rows × 3 columns

	Y	Year	Corp
	Int32	Int32	Cat…
1	0	75	G
2	2	76	G
3	2	77	G
4	1	78	G
5	0	79	G
6	0	80	G
7	1	81	G
8	1	82	G
9	0	83	G
10	3	84	G
11	0	85	G
12	2	86	G
13	1	87	G
14	0	88	G
15	0	89	G
16	1	90	G
17	0	91	G
18	1	92	G
19	0	93	G
20	1	94	G
21	0	75	I
22	0	76	I
23	0	77	I
24	2	78	I
25	0	79	I
26	3	80	I
27	0	81	I
28	2	82	I
29	0	83	I
30	0	84	I
⋮	⋮	⋮	⋮

# Año
x = sort(unique(prussian[!,:Year]))

20-element Array{Int32,1}:
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94

# Número de muertes por patada de caballo
y = [ sum(prussian[!,:Y][ prussian[!,:Year] .== x[i] ]) for i in 1:length(x) ]

20-element Array{Int64,1}:
  3
  5
  7
  9
 10
 18
  6
 14
 11
  9
  5
 11
 15
  6
 11
 17
 12
 15
  8
  4

plot(p)

plot(p)

Simulación de la cola de una distribución

Sea $s\in\mathbb{R}$ y $X\sim \mathcal{L}\left( f(x) \right)$ con $f$ una función de densidad de probabilidad en $\mathbb{R}$ que es monótona en $[s,\infty)$.

Consideremos $U$ una variable aleatoria que toma valores en $(0,1]$, $\sigma >0$ y $Y=g(U)=s-\sigma\log(U)$ entonces $g^{-1}(y)=e^{ (s-y)/\sigma}$ y $U$ tiene densidad

\begin{align*} f_U(u) & = \frac{f_Y(y)}{\left| \frac{\mathrm{d}u}{\mathrm{d}y} \right|} = \frac{f_Y\left( y\right)}{\left| \frac{\mathrm{d}g^{-1}(y)}{\mathrm{d}y} \right|} = \frac{f_Y(y)}{\left| \frac{1}{g'\left( g^{-1}(y) \right) } \right|} = f_Y(y)\left| g'\left( g^{-1}(y) \right) \right| \\ & = f_Y(y) \left| -\frac{\sigma}{g^{-1}(y)} \right| =\sigma f_Y(y)e^{ (y-s)/ \sigma} = \frac{ \sigma f_Y(y) }{ u } \end{align*}

Por construcción la variable aleatoria $Y$ toma valores en $[s,\infty )$. Por lo que podemos determinar la ley de probabilidad de $U$ al elegir $f_Y(y)=f$, nuestra densidad de probabilidad de interés. Así las cosas, si $(Z_1,Z_2)$ es tomado uniformemente debajo de la gráfica de $f_U(u)$ entonces $Z_1$ sería una realización de $\mathcal{L}(f)$ restringida a $[s,\infty )$, es decir usamos aceptación-rechazo. Este método es especialmente conveniente si con $f_Y=f$ se tiene que $f_U$ es estrictamente decreciente en $[s,\infty)$ ya que entonces podemos usar la cota $Mg(x)=f_U(s)$.

Ejemplo

Si $f(x)=\frac{1}{\sqrt{2\pi}}e^{-x^2/2}$ entonces $f_U(u)=\frac{ \sigma }{\sqrt{2\pi}}e^{-y^2/2}e^{(y-s)/\sigma} $ es estrictamente decreciente si y sólo si

\begin{align*} \frac{\mathrm{d}}{\mathrm{d}y}\left( e^{-y^2/2 + (y-s)/\sigma} \right) <0 & \iff e^{-y^2/2 + (y-s)/\sigma}\left( -y + 1/\sigma \right) < 0 \\ & \iff \frac{1}{y} < \sigma. \end{align*}

Recordemos que $s < y$ por lo que se puede elegir $\sigma\in [1/s, \infty )$. Así las cosas podemos utilizar $Mg(u)=\sigma e^{-s^2/2}/\sqrt{2\pi}$ y la regla de aceptación

\begin{align*} U_2 \leq \frac{f_U(U_1)}{f_U(s)}=\frac{ e^{-Y^2/2 + s^2/2} }{U_1}=e^{-Y^2/2 + s^2/2 + (Y-s)/\sigma}. \end{align*}

La probabilidad de aceptación en este caso es un factor de $1/\sigma$ al ser $\sigma$ factor de $Mg$. Por esta razón tomamos el mínimo valor posible de $\sigma$ que es $1/s$.

Simulación de una Normal restringida a la cola $[s,\infty )$

1) Genera $U_1, U_2 \sim \text{Uniforme}(0,1)$ independientes y $Y=s - \frac{1}{s}\log(U_1)$.

2) Si $$U_2 \leq e^{-Y^2/2 + s^2/2 +s(Y-s) }$$ aceptamos $X=Y$, en caso contrario regresamos a 1).

Observando que $$U_2 \leq e^{-Y^2/2 + s^2/2 +s(Y-s) } \iff -2 log(U_2) \leq \left( \frac{\log(U_1)}{s} \right)^2$$ obtenemos un algoritmo menos costoso:

1) Genera $U_1, U_2 \sim \text{Uniforme}(0,1)$ independientes y $Y= -\frac{1}{s}\log(U_1)$.

2) Si $$-2\log(U_2) > Y^2$$ aceptamos $X=s+Y$, en caso contrario regresamos a 1).